home *** CD-ROM | disk | FTP | other *** search
/ L' Effet Pommier 3 / L'Effet Pommier - Volume 03.iso / Programmation / gray image 2.1 / image.cc < prev    next >
Text File  |  1995-10-08  |  22KB  |  794 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *               Grayscale Image
  6.  *         Implementation of the Primitive Operations
  7.  *
  8.  *   The image is represented as a Pixmap, i.e. a matrix of pixels
  9.  *    each of them specifies the gray level at a particular point
  10.  *
  11.  ************************************************************************
  12.  */
  13.  
  14. #pragma implementation "image.h"
  15. #include "image.h"
  16. #include "std.h"
  17. #include <minmax.h>
  18. #include <math.h>
  19.  
  20.  
  21.  
  22. /*
  23.  *------------------------------------------------------------------------
  24.  *            Constructors and destructors
  25.  */
  26.  
  27. void IMAGE::allocate(
  28.     const card no_rows,        // No. of rows
  29.     const card no_cols,        // No. of cols
  30.     const card depth        // No. of bits per pixel
  31. )
  32. {
  33.   valid_code = IMAGE_val_code;
  34.  
  35.   assure((ncols=no_cols) > 0, "Zero image width unexpected");
  36.   assure((nrows=no_rows) > 0, "Zero image height unexpected");
  37.   assure((bits_per_pixel=depth) > 0 && depth <= GRAY_MAXBIT, 
  38.      "Zero or too large no. of bits per pixel");
  39.  
  40.   name    = "";
  41.   npixels = nrows * ncols;
  42.  
  43.   assert( (scanrows = (GRAY **)calloc(nrows,sizeof(GRAY *))) != 0 );
  44.   assert( (pixels   = (GRAY *)calloc(npixels,sizeof(GRAY))) != 0 );
  45.  
  46.   register GRAY * pp = pixels;
  47.   register GRAY **sp = scanrows;
  48.   while( pp < pixels + npixels )    // Build the row index
  49.     *sp++ = pp, pp += ncols;
  50.   assert( sp == scanrows + nrows );
  51. }
  52.  
  53.                 // Set a new image name
  54. void IMAGE::set_name(const char * new_name)
  55. {
  56.   if( name != 0 && name[0] != '\0' )    // Dispose of the previous image name
  57.     delete name;
  58.  
  59.   if( new_name == 0 || new_name[0] == '\0' )
  60.     name = "";                // Image is anonymous now
  61.   else
  62.     name = new char[strlen(new_name)+1], strcpy(name,new_name);
  63. }
  64.  
  65.                     // Routing constructor module
  66. IMAGE::IMAGE(const IMAGE_CREATORS_1op op, const IMAGE& prototype)
  67. {
  68.   switch(op)
  69.   {
  70.     case Expand:
  71.          _expand(prototype);
  72.      break;
  73.  
  74.     case Shrink:
  75.      _shrink(prototype);
  76.      break;
  77.  
  78.     default:
  79.      _error("Operation %d is not yet implemented",op);
  80.   }
  81. }
  82.  
  83.  
  84. IMAGE::~IMAGE(void)        // Dispose the image struct
  85. {
  86.   is_valid();
  87.   if( name != 0 && name[0] != '\0' )
  88.     delete name;
  89.  
  90.   free(scanrows);
  91.   free(pixels);
  92.   valid_code = 0;
  93. }
  94.  
  95.  
  96. /*
  97.  *------------------------------------------------------------------------
  98.  *             Single Image operations
  99.  */
  100.  
  101.                 // Clip pixel values to
  102.                 // [0,1<<bits_per_pixel-1].
  103.                 // A pixel with the value outside that range
  104.                 // (i.e. negative or too big) is set to
  105.                 // 0 or 1<<bits_per_pixel-1, resp.
  106. IMAGE& IMAGE::clip_to_intensity_range(void)
  107. {
  108.   is_valid();
  109.   const int maxval = (1<<bits_per_pixel)-1;
  110.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  111.   for(; cp < (GRAY_SIGNED *)pixels + npixels; cp++)
  112.     if( *cp < 0 )
  113.       *cp = 0;
  114.     else if( *cp > maxval )
  115.       *cp = maxval;
  116.  
  117.   return *this;
  118. }
  119.  
  120.                 // Perform a transformation
  121.                 // pixel = |(signed)pixel|
  122.                 // on all the pixels of the image
  123. IMAGE& IMAGE::abs(void)
  124. {
  125.   is_valid();
  126.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  127.   for(; cp < (GRAY_SIGNED *)pixels + npixels; cp++)
  128.     if( *cp < 0 )
  129.       *cp = -(*cp);
  130.  
  131.   return *this;
  132. }
  133.  
  134.                 // Apply a user-defined action to each pixel
  135.                 // while image is efficiently traversed
  136.                 // row-by-row
  137. IMAGE& IMAGE::apply(PixelAction& action)
  138. {
  139.   is_valid();
  140.   action.nrows = nrows, action.ncols = ncols;
  141.   register GRAY * cp = pixels;
  142.   for(action.row=0; action.row<nrows; action.row++)
  143.     for(action.col=0; action.col<ncols; action.col++)
  144.       action.operation(*cp++);
  145.  
  146.   return *this;
  147. }
  148.  
  149.                 // Perform the histogram equalization
  150. IMAGE& IMAGE::equalize(const int no_grays)
  151. {
  152.   is_valid();
  153.   int orig_no_grays = 1 << bits_per_pixel;
  154.   assert( no_grays <= orig_no_grays );
  155.  
  156.   int * const histogram = (int *)calloc(orig_no_grays,sizeof(int));
  157.  
  158.                 // Evaluate the histogram of an image
  159.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  160.   for(; cp < (GRAY_SIGNED *)pixels + npixels; cp++)
  161.   {
  162.     if( *cp <= 0 )
  163.       *cp = 0;
  164.     else if( *cp >= orig_no_grays )
  165.       *cp = orig_no_grays-1;
  166.     histogram[*cp]++;
  167.   }
  168.  
  169.   const int pixels_per_bin_optimal = (npixels + no_grays - 1)/no_grays;
  170.   const int gray_shade_subsample_factor = 
  171.     (orig_no_grays + no_grays -1)/no_grays;
  172.   register int pixel;
  173.   int new_pixel, new_pixel_old;
  174.   int accumulated = 0;
  175.   int optimal_accumulated = pixels_per_bin_optimal;
  176.                     // Mapping from pixels of original
  177.                       // image to [new_pixel_old,new_pixel]
  178.   short * const look_up_table = (short *)malloc(orig_no_grays*sizeof(short));
  179.                     // (actually, just to the center of
  180.                     // this interval)
  181.  
  182.                       // Equalizing the histogram
  183.   for(pixel=0,new_pixel=0,new_pixel_old=0; pixel<orig_no_grays; pixel++)
  184.   {
  185.     accumulated += histogram[pixel];
  186.     while( accumulated > optimal_accumulated )
  187.       new_pixel += gray_shade_subsample_factor, 
  188.       optimal_accumulated += pixels_per_bin_optimal;
  189.     assert( new_pixel < orig_no_grays );
  190.     look_up_table[pixel] = (new_pixel > 
  191.                 new_pixel_old+gray_shade_subsample_factor ?
  192.                 (new_pixel + new_pixel_old)/2 :
  193.                 new_pixel);
  194.     new_pixel_old = new_pixel;
  195.   }
  196.  
  197.                 // Update the image according to the LUT
  198.   for(cp = (GRAY_SIGNED *)pixels; cp < (GRAY_SIGNED *)pixels + npixels; cp++)
  199.     *cp = look_up_table[*cp];
  200.  
  201.   free(histogram);
  202.   free(look_up_table);
  203.   
  204.   return *this;
  205. }
  206.  
  207.                 // Compute the 1. norm of the entire image
  208.                 // SUM{ |(signed)pixel[i,j]| }
  209. double IMAGE::norm_1(void) const
  210. {
  211.   is_valid();
  212.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  213.   double sum = 0;
  214.   for(; cp < (GRAY_SIGNED *)pixels + npixels; )
  215.     sum += ::abs(*cp++);
  216.  
  217.   return sum;
  218. }
  219.  
  220.                 // Compute the square of the 2. norm of 
  221.                 // the entire image
  222.                 // SUM{ |(signed)pixel[i,j]|^2 }
  223. double IMAGE::norm_2_sqr(void) const
  224. {
  225.   is_valid();
  226.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  227.   register double sum = 0;
  228.   for(; cp < (GRAY_SIGNED *)pixels + npixels; )
  229.     sum += ::sqr((long int)*cp++);
  230.  
  231.   return sum;
  232. }
  233.  
  234.                 // Compute the infinity norm of the
  235.                 // entire image
  236.                 // MAX{ |(signed)pixel[i,j]| }
  237. int IMAGE::norm_inf(void) const
  238. {
  239.   is_valid();
  240.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  241.   register int maxp = 0;
  242.   for(; cp < (GRAY_SIGNED *)pixels + npixels; )
  243.     maxp = ::max(::abs(*cp++),maxp);
  244.  
  245.   return maxp;
  246. }
  247.  
  248.                 // Find extremum values of image pixels
  249. Extrema::Extrema(const IMAGE& image)
  250.     : max_pixel(0,0), min_pixel(0,0)
  251. {
  252.   image.is_valid();
  253.   max_value = min_value = image(0,0);
  254.  
  255.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)image.pixels;
  256.   register int i,j;
  257.   for(i=0; i<image.nrows; i++)
  258.     for(j=0; j<image.ncols; j++, cp++)
  259.       if( *cp > max_value )
  260.     max_value = *cp, max_pixel.row_val = i, max_pixel.col_val = j;
  261.       else if( *cp < min_value )
  262.     min_value = *cp, min_pixel.row_val = i, min_pixel.col_val = j;
  263. }
  264.  
  265.                 // Normalize pixel values to be
  266.                 // in range 0..1<<bits_per_pixel-1
  267. IMAGE& IMAGE::normalize_for_display(void)
  268. {
  269.   is_valid();
  270.   Extrema extrema(*this);
  271.   message("\nImage_normalization:"
  272.       "\n\tmin pixel value %d, max pixel value %d",
  273.       extrema.min(), extrema.max());
  274.   if( extrema.max() != extrema.min() )
  275.   {
  276.     double factor = (double)((1<<bits_per_pixel)-1) /
  277.             ( extrema.max() - extrema.min() );
  278.     message("\n\tNormalization is as follows: (pixel - %d)*%g\n",
  279.         extrema.min(), factor);
  280.     register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  281.     for(; cp < (GRAY_SIGNED *)pixels + npixels; cp++)
  282.       *cp = (int)( (*cp - extrema.min()) * factor + 0.5 );
  283.   }
  284.  
  285.   return *this;
  286. }
  287.  
  288.                 // Print some info about the image
  289. void IMAGE::info(void) const
  290. {
  291.   message("\nimage %dx%dx%d '%s' ",nrows,ncols,bits_per_pixel,name);
  292. }
  293.                 // Print the image as a table of pixel values
  294.                 // (zeros are printed as dots)
  295. void IMAGE::print(const char * title) const
  296. {
  297.   is_valid();
  298.   message("\nImage %dx%dx%d '%s' is as follows\n\n",nrows,ncols,
  299.       bits_per_pixel,title);
  300.  
  301.   register int i,j;
  302.   for(i=0; i<nrows; i++)
  303.   {
  304.     for(j=0; j<ncols; j++)
  305.       if( (*this)(i,j) == 0 )
  306.     message("   . ");
  307.       else
  308.     message("%4d ",(GRAY_SIGNED)(*this)(i,j));
  309.     message("\n");
  310.   }
  311.   message("Done\n");
  312. }
  313.  
  314. /*
  315.  *------------------------------------------------------------------------
  316.  *             Image-scalar operations
  317.  *       Check to see if the preedicate "(signed)pixel operation scalar"
  318.  *        holds for ALL the pixels of the image
  319.  */
  320.  
  321.                 // is "(signed)pixel OP val" true for all
  322.                 // pixels?
  323.  
  324. #define COMPARISON_WITH_SCALAR(OP)                    \
  325.                                     \
  326. bool IMAGE::operator OP (const int val) const                \
  327. {                                    \
  328.   is_valid();                                \
  329.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;            \
  330.   for(; cp < (GRAY_SIGNED *)pixels + npixels; )                \
  331.     if( !(*cp++ OP val) )                        \
  332.       return false;                            \
  333.                                     \
  334.   return true;                                \
  335. }                                    \
  336.  
  337. COMPARISON_WITH_SCALAR(==)
  338. COMPARISON_WITH_SCALAR(!=)
  339. COMPARISON_WITH_SCALAR(<)
  340. COMPARISON_WITH_SCALAR(<=)
  341. COMPARISON_WITH_SCALAR(>)
  342. COMPARISON_WITH_SCALAR(>=)
  343.  
  344. #undef COMPARISON_WITH_SCALAR
  345.  
  346.  
  347. /*
  348.  *------------------------------------------------------------------------
  349.  *             Image-scalar operations
  350.  *          Modify every pixel according to the operation
  351.  */
  352.  
  353.                 // For every pixel, do `pixel OP value`
  354. #define COMPUTED_VAL_ASSIGNMENT(OP)                    \
  355.                                     \
  356. IMAGE& IMAGE::operator OP (const int val)                \
  357. {                                    \
  358.   is_valid();                                \
  359.   register GRAY * cp = pixels;                        \
  360.   while( cp < pixels+npixels )                        \
  361.     *cp++ OP val;                            \
  362.                                     \
  363.   return *this;                                \
  364. }                                    \
  365.  
  366. COMPUTED_VAL_ASSIGNMENT(=)
  367. COMPUTED_VAL_ASSIGNMENT(+=)
  368. COMPUTED_VAL_ASSIGNMENT(-=)
  369. COMPUTED_VAL_ASSIGNMENT(*=)
  370. COMPUTED_VAL_ASSIGNMENT(|=)
  371. COMPUTED_VAL_ASSIGNMENT(&=)
  372. COMPUTED_VAL_ASSIGNMENT(^=)
  373.  
  374. #undef COMPUTED_VAL_ASSIGNMENT
  375.  
  376.                 // Shift the value of all the pixels
  377. IMAGE& IMAGE::operator <<= (const int val)
  378. {
  379.   is_valid();
  380.   if( ::abs(val) >= GRAY_MAXBIT )
  381.     _error("Very fishy shift factor: %d",val);
  382.  
  383.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  384.   while( cp < (GRAY_SIGNED *)pixels+npixels )
  385.     *cp++ <<= val;
  386.  
  387.   return *this;
  388. }
  389.  
  390.                 // Shift the value of all the pixels
  391. IMAGE& IMAGE::operator >>= (const int val)
  392. {
  393.   is_valid();
  394.   if( ::abs(val) >= GRAY_MAXBIT )
  395.     _error("Very fishy shift factor: %d",val);
  396.  
  397.   register GRAY_SIGNED * cp = (GRAY_SIGNED *)pixels;
  398.   while( cp < (GRAY_SIGNED *)pixels+npixels )
  399.     *cp++ >>= val;
  400.  
  401.   return *this;
  402. }
  403.  
  404. /*
  405.  *------------------------------------------------------------------------
  406.  *             Image-Image operations
  407.  *          Modify the target image according to the operation
  408.  */
  409.  
  410. IMAGE& IMAGE::operator = (const IMAGE& source)
  411. {
  412.   are_compatible(*this,source);
  413.   memcpy(pixels,source.pixels,npixels*sizeof(GRAY));
  414.   return *this;
  415. }
  416.  
  417. bool operator == (const IMAGE& im1, const IMAGE& im2)
  418. {
  419.   are_compatible(im1,im2);
  420.   return (memcmp(im1.pixels,im2.pixels,im1.npixels*sizeof(GRAY)) == 0);
  421. }
  422.  
  423.                 // Do "image OP src"
  424. #define COMPUTED_ASSIGNMENT(OP)                        \
  425.                                     \
  426. IMAGE& IMAGE::operator OP (const IMAGE& src)                \
  427. {                                    \
  428.   are_compatible(*this,src);                        \
  429.                                     \
  430.   register GRAY * sp = src.pixels;                    \
  431.   register GRAY * tp = pixels;                        \
  432.   while( tp < pixels + npixels )                    \
  433.     *tp++ OP *sp++;                            \
  434.                                     \
  435.   return *this;                                \
  436. }                                    \
  437.  
  438. COMPUTED_ASSIGNMENT(+=)
  439. COMPUTED_ASSIGNMENT(-=)
  440. COMPUTED_ASSIGNMENT(|=)
  441. COMPUTED_ASSIGNMENT(&=)
  442. COMPUTED_ASSIGNMENT(^=)
  443.  
  444. #undef COMPUTED_ASSIGNMENT
  445.  
  446. #if 0
  447.                 // Modified addition
  448.                 //    Target += scalar*Source
  449. IMAGE& add(IMAGE& target, const int scalar,const IMAGE& source)
  450. {
  451.   are_compatible(target,source);
  452.  
  453.   register GRAY * sp = source.pixels;
  454.   register GRAY * tp = target.pixels;
  455.   while( tp < target.pixels + target.npixels )
  456.     *tp++ += scalar * *sp++;
  457.   
  458.   return target;
  459. }
  460.                 // Shift the source to 'pos', clip it 
  461.                 // if necessary, multiply by the scalar,
  462.                 // and add
  463. IMAGE& IMAGE::shift_clip_add(rowcol pos, const int scalar, const IMAGE& source)
  464. {
  465.   source.is_valid();
  466.   is_valid();
  467.  
  468.                     // Find an intersection of 'source'
  469.                     // shifted to 'pos' (relative to
  470.                     // the image) with the (target) image
  471.                     // Only this rectangle will be 
  472.                     // affected by addition
  473.   rowcol s_orig(max(0,-pos.row()),max(0,-pos.col()));
  474.   rowcol t_orig(max(0,pos.row()),max(0,pos.col()));
  475.                     // No. of rows in the intersection
  476.   const int irows = min(nrows - t_orig.row(),source.nrows - s_orig.row());
  477.                     // No. of cols in the intersection
  478.   const int icols = min(ncols - t_orig.col(),source.ncols - s_orig.col());
  479.   if( irows <= 0 || icols <= 0 )
  480.     source.info(),
  481.     info(),
  482.     _error("The first image shifted by (%d,%d) does not intersect the other",
  483.        pos.row(),pos.col());
  484.  
  485.                     // See the explanation in the class
  486.                     // Rectangle
  487.   const int s_inc_to_nextrow = source.ncols - icols;    
  488.   const int t_inc_to_nextrow = ncols - icols;    
  489.  
  490.   register GRAY_SIGNED * sp = 
  491.     (GRAY_SIGNED *)&(source.scanrows[s_orig.row()])[s_orig.col()];
  492.   register GRAY_SIGNED * tp = 
  493.     (GRAY_SIGNED *)&(scanrows[t_orig.row()])[t_orig.col()];
  494.  
  495.   register int i,j;
  496.   for(i=0; i<irows; i++, sp += s_inc_to_nextrow, tp += t_inc_to_nextrow)
  497.     for(j=0; j<icols; j++)
  498.       *tp++ += scalar * *sp++;
  499.  
  500.   return *this;
  501. }
  502. #endif
  503.  
  504.                 // Evaluate the scalar product of two images
  505. double operator * (const IMAGE& im1, const IMAGE& im2)
  506. {
  507.   are_compatible(im1,im2);
  508.  
  509.   register GRAY_SIGNED * p1 = (GRAY_SIGNED *)im1.pixels;
  510.   register GRAY_SIGNED * p2 = (GRAY_SIGNED *)im2.pixels;
  511.   register double sum = 0;
  512.  
  513.   while( p1 < (GRAY_SIGNED *)im1.pixels + im1.npixels )
  514.     sum += *p1++ * *p2++;
  515.  
  516.   return sum;
  517. }
  518.  
  519.             // Estimate the norm of the difference 
  520.             // between two (signed) image
  521.                 // SUM{ |(signed)pixel[i,j]| }
  522. double norm_1(const IMAGE& im1, const IMAGE& im2)
  523. {
  524.   are_compatible(im1,im2);
  525.  
  526.   register GRAY_SIGNED * p1 = (GRAY_SIGNED *)im1.pixels;
  527.   register GRAY_SIGNED * p2 = (GRAY_SIGNED *)im2.pixels;
  528.   register double sum = 0;
  529.  
  530.   while( p1 < (GRAY_SIGNED *)im1.pixels + im1.npixels )
  531.     sum += ::abs(*p1++ - *p2++);
  532.  
  533.   return sum;
  534. }
  535.  
  536.                 // SUM{ |(signed)pixel[i,j]|^2 }
  537. double norm_2_sqr(const IMAGE& im1, const IMAGE& im2)
  538. {
  539.   are_compatible(im1,im2);
  540.  
  541.   register GRAY_SIGNED * p1 = (GRAY_SIGNED *)im1.pixels;
  542.   register GRAY_SIGNED * p2 = (GRAY_SIGNED *)im2.pixels;
  543.   register double sum = 0;
  544.  
  545.   while( p1 < (GRAY_SIGNED *)im1.pixels + im1.npixels )
  546.     sum += ::sqr((long int)(*p1++ - *p2++));
  547.  
  548.   return sum;
  549. }
  550.                 // MAX{ |(signed)pixel[i,j]| }
  551. int norm_inf(const IMAGE& im1, const IMAGE& im2)
  552. {
  553.   are_compatible(im1,im2);
  554.  
  555.   register GRAY_SIGNED * p1 = (GRAY_SIGNED *)im1.pixels;
  556.   register GRAY_SIGNED * p2 = (GRAY_SIGNED *)im2.pixels;
  557.   register int maxp = 0;
  558.  
  559.   while( p1 < (GRAY_SIGNED *)im1.pixels + im1.npixels )
  560.     maxp = max(maxp,::abs(*p1++ - *p2++));
  561.  
  562.   return maxp;
  563. }
  564.  
  565.  
  566. void compare(            // Compare the two images
  567.     const IMAGE& image1,    // and print out the result of comparison
  568.     const IMAGE& image2,
  569.     const char * title )
  570. {
  571.   register int i,j;
  572.  
  573.   are_compatible(image1,image2);
  574.  
  575.   message("\n\nComparison of two images %dx%dx%d\n%s\n",image1.nrows,
  576.      image1.ncols,image1.bits_per_pixel,title);
  577.  
  578.   if( image1.bits_per_pixel != image2.bits_per_pixel )
  579.     message("\nNote, images have different depth, %d and %d\n",
  580.          image1.bits_per_pixel, image2.bits_per_pixel);
  581.  
  582.   double norm1 = 0, norm2 = 0;        // Norm of the images
  583.   double ndiff = 0;            // Norm of the difference
  584.   int imax=0,jmax=0,difmax = -1;    // For the pixels that differ most
  585.                     // Image scanline pointer
  586.   register GRAY_SIGNED *rowp1 = (GRAY_SIGNED *)image1.pixels;
  587.   register GRAY_SIGNED *rowp2 = (GRAY_SIGNED *)image2.pixels;
  588.  
  589.   bool im1_neg_pixel = false;        // Flags whether negative pixels
  590.   bool im2_neg_pixel = false;        // were encountered
  591.  
  592.   for(i=0; i < image1.nrows; i++)
  593.     for(j=0; j < image1.ncols; j++)
  594.     {
  595.       int pv1 = *rowp1++;
  596.       int pv2 = *rowp2++;
  597.       int diff = abs(pv1-pv2);
  598.  
  599.       if( pv1 < 0 )
  600.     im1_neg_pixel = true, pv1 = -pv1;
  601.       if( pv2 < 0 )
  602.     im2_neg_pixel = true, pv2 = -pv2;
  603.  
  604.       if( diff > difmax )
  605.       {
  606.     difmax = diff;
  607.     imax = i;
  608.     jmax = j;
  609.       }
  610.       norm1 += pv1;
  611.       norm2 += pv2;
  612.       ndiff += diff;
  613.     }
  614.  
  615.   if( im1_neg_pixel )
  616.     message("\n*** Warning: a pixel with negative value has been encountered "
  617.         "in image 1\n");
  618.   if( im2_neg_pixel )
  619.     message("\n*** Warning: a pixel with negative value has been encountered "
  620.         "in image 2\n");
  621.  
  622.   message("\nMaximal discrepancy    \t\t%d",difmax);
  623.   message("\n   occured at the pixel\t\t(%d,%d)",imax,jmax);
  624.   const int pv1 = image1(imax,jmax);
  625.   const int pv2 = image2(imax,jmax);
  626.   message("\n Image 1 pixel is    \t\t%d",pv1);
  627.   message("\n Image 2 pixel is    \t\t%d",pv2);
  628.   message("\n Absolute error v2[i]-v1[i]\t\t%d",pv2-pv1);
  629.   message("\n Relative error\t\t\t\t%g\n",
  630.      (pv2-pv1)/max((float)abs(pv2+pv1)/2,1e-7) );
  631.  
  632.   message("\nL1 norm of image 1 per pixel \t\t%g",
  633.                       (double)norm1/image1.npixels);
  634.   message("\nL1 norm of image 2 per pixel \t\t%g",
  635.                     (double)norm2/image2.npixels);
  636.   message("\nL1 norm of image1-image2 per pixel\t\t\t%g",
  637.                     (double)ndiff/image1.npixels);
  638.   message("\n||Image1-Image2||/sqrt(||Image1|| ||Image2||)\t%g\n\n",
  639.       ndiff/max( sqrt(norm1*norm2), 1e-7 )         );
  640.  
  641. }
  642.  
  643. /*
  644.  *------------------------------------------------------------------------
  645.  *                Expansion/Shrinking of the image
  646.  */
  647.  
  648.                 // Expand the prototype twice in each
  649.                 // dimension
  650.                 // In other words, blow each pixel of the
  651.                 // prototype to the 2x2 square
  652. void IMAGE::_expand(const IMAGE& prototype)
  653. {
  654.   prototype.is_valid();
  655.   allocate(2*prototype.nrows,2*prototype.ncols,prototype.bits_per_pixel);
  656.  
  657.   register GRAY * tp = pixels;
  658.   register GRAY * pp = prototype.pixels;
  659.   register int i,j;
  660.  
  661.   for(i=0; i<prototype.nrows; i++, tp+=ncols)
  662.     for(j=0; j<prototype.ncols; j++)
  663.     {
  664.       register int pixel = *pp++;
  665.       *(tp+ncols) = pixel;        // Fill out 2x2 square of the
  666.       *tp++       = pixel;        // blown image with a pixel value
  667.       *(tp+ncols) = pixel;
  668.       *tp++       = pixel;
  669.     }
  670.  
  671.   assert( pp == prototype.pixels + prototype.npixels );
  672.   assert( tp == pixels + npixels );
  673. }
  674.  
  675.                 // Shrink the prototype twice in each
  676.                 // dimension. The image row and column
  677.                 // dimensions are assumed to be even numbers
  678.                 // Image is shrunk by replacing 2x2 square
  679.                 // by one pixel with an average intensity
  680.                 // over the square
  681. void IMAGE::_shrink(const IMAGE& prototype)
  682. {
  683.   prototype.is_valid();
  684.   if( (prototype.nrows & 1) || (prototype.ncols & 1) )
  685.     _error("No of rows and columns in the image to shrink should both be even",
  686.        (prototype.info(),0));
  687.   allocate(prototype.nrows/2,prototype.ncols/2,prototype.bits_per_pixel);
  688.  
  689.   register GRAY * tp = pixels;
  690.   register GRAY * pp = prototype.pixels;
  691.   register int i,j;
  692.  
  693.   for(i=0; i<nrows; i++)
  694.   {                    // Average row-by-row
  695.     for(j=0; j<ncols; j++)
  696.       *tp = *pp++, *tp++ += *pp++;    // Sum of pairs of the prototype pixels
  697.     tp -= ncols;
  698.     for(j=0; j<ncols; j++, tp++, pp+=2)
  699.     {
  700.       register int av = *tp + pp[0] + pp[1];    // This is the aver of 4 pixels
  701.       *tp = ( av & 2 ? (av >> 2) + 1 : av >> 2 );    // /4 with rounding
  702.     }
  703.   }
  704.   assert( pp == prototype.pixels + prototype.npixels );
  705.   assert( tp == pixels + npixels );
  706. }
  707.  
  708.  
  709.                 // Assign another to '*this' resizing
  710.                 // '*this' as necessary to fit exactly
  711.                 // Squeezing is done by subsampling (dropping
  712.                 // some image pixels), while stretching is done
  713.                 // by repeating pixels/scanlines
  714.                 // Note, the ratio of sizes of the images
  715.                 // can be *anything*
  716. IMAGE& IMAGE::coerce(const IMAGE& another)
  717. {
  718.                     // Make the width even
  719. //  width  = ( dest_rect.q_width() + 1 ) & (~1);
  720. //  height = dest_rect.q_height();
  721. //  assert( width > 0 && height > 0 );
  722.  
  723.                     // Figure out how much to
  724.                     // sqeeze/stretch
  725.   const float ratio_h = (float)ncols/another.ncols;
  726.   const float ratio_v = (float)nrows/another.nrows;
  727.  
  728.   float curr_v = 0, curr_h;        // Precise curr. loc. in window coord
  729.   card dest_y = 0, dest_x;        // approx (rounded-off) coordinates
  730.   register GRAY * drp = pixels;        // Scanline ptr for fitted image
  731.   register const GRAY * spp = another.pixels; // Scanline of 'another' image
  732.  
  733.                      // Note, dest_x,dest_y is incremented
  734.                     // by (inc_v,inc_h), but curr_v is
  735.                     // incremented by ratio_v, etc.
  736.                     // Draw horizontally (row-by-row)
  737.   while( spp < another.pixels + another.npixels )
  738.   {
  739.     curr_v += ratio_v;
  740.     const GRAY * const spp_end = spp + another.ncols; // end of scanline
  741.  
  742.     int inc_v = (int)curr_v - dest_y;
  743.     if( inc_v == 0 )            // inc_v = 0 means this row of 
  744.     {                    // orig_image is to be skipped
  745.       spp = spp_end;
  746.       continue;
  747.     }
  748.     dest_y = (int)curr_v;
  749.  
  750.     curr_h = 0;
  751.     dest_x = 0;
  752.     register GRAY * dpp = drp;        // Squeeze/stretch pixels in a scanline
  753.     for(; spp < spp_end; spp++)
  754.     {
  755.       curr_h += ratio_h;
  756.  
  757.       int inc_h = (int)curr_h - dest_x;
  758.       dest_x = (int)curr_h;
  759.       while( inc_h-- > 0 )
  760.     *dpp++ = *spp;
  761.     }
  762.     if( dpp == drp+ncols-1 )        // Can happen due to roundoff arithm
  763.       *dpp++ = *(spp-1);
  764.     assert( dpp == drp+ncols );
  765.                     // Stretching vertically involves
  766.                     // duplicating scanlines
  767.     while( --inc_v > 0 )
  768.       memcpy((void *)(drp+ncols),(void *)drp,ncols*sizeof(GRAY)),
  769.       drp += ncols;
  770.     drp += ncols;
  771.   }
  772.   if( drp == pixels + npixels-ncols )    // Duplicate the last scanline
  773.     memcpy((void *)drp,(void *)(drp-ncols),ncols*sizeof(GRAY)), // in case
  774.     drp += ncols;            // of roundoff snag
  775.  
  776.   assert( drp == pixels + npixels );
  777.  
  778.   return *this;
  779. }
  780.  
  781.  
  782.                     // Verify the pixels have the value
  783.                     // that is expected
  784. void verify_pixel_value(const IMAGE& im, const GRAY val)
  785. {
  786.   register int i,j;
  787.   for(i=0; i<im.q_nrows(); i++)
  788.     for(j=0; j<im.q_ncols(); j++)
  789.       if( im(i,j) != val )
  790.     _error("Pixel [%d,%d] has the value 0x%x different from expected 0x%x",
  791.            i,j,im(i,j),val);
  792. }
  793.  
  794.